
/*
 * Implementamos las funciones definidas en inversa.h
 * Actualmente la única función definida es "inversa".
 */

#include "inversa.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>



/* Definición de tipos de datos */

/*
 * Esta estructura la utilizamos para indicar cual fue el resultado de la
 * factorzación LU.
 */
struct RESULTADO_FACTORIZACION {
    struct MATRIZ *matrizL;
    struct MATRIZ *matrizU;
};



/* Prototipado de funciones locales */
int analizarTriangularizacion(struct MATRIZ *matriz);
struct RESULTADO_FACTORIZACION triangularizar(struct MATRIZ *matriz);
struct MATRIZ *invertirMatrizL(struct MATRIZ *matrizL);
struct MATRIZ *invertirMatrizU(struct MATRIZ *matrizU);
int sustitucionDirecta(struct MATRIZ *coefs, struct VECTOR *ti, struct VECTOR *solucion);
int sustitucionInversa(struct MATRIZ *coefs, struct VECTOR *ti, struct VECTOR *solucion);



/* Implementación de las funciones exportadas y locales */

/*
 * Implementación de la función definida en cabecera.
 */
struct MATRIZ* inversa(struct MATRIZ *matriz) {

    // Si no se puede triangularizar, devolvemos NULL
    if(!analizarTriangularizacion(matriz)) {
        printf("\nERROR: La matriz no se puede triangularizar");
        return NULL;
    }

    // Realizamos la factorizacion LU de la matriz
    struct RESULTADO_FACTORIZACION resultadoLU = triangularizar(matriz);

    // Invertimos la matriz L y eliminamos L
    struct MATRIZ *inversaL = invertirMatrizL(resultadoLU.matrizL);
    eliminarMatriz(resultadoLU.matrizL);

    // Invertimos la matriz U y eliminamos U
    struct MATRIZ *inversaU = invertirMatrizU(resultadoLU.matrizU);
    eliminarMatriz(resultadoLU.matrizU);

    // Si algunas de las matrices es nula, devolver NULL
    if(inversaL == NULL || inversaU == NULL) {
        printf("\nNo se pudo calcular la inversa de la matriz");
        return NULL;
    }

    // Multiplicamos ambas inversas y obtenemos la inversa de la matriz
    struct MATRIZ *matrizInvertida = multiplicarMatricesLU(inversaU,inversaL);
    
    // Eliminamos las inversas de L y U
    eliminarMatriz(inversaL);
    eliminarMatriz(inversaU);

    // Devolvemos la inversa de la matriz
    return matrizInvertida;

}


/*
 * Verificamos que no existan ceros en la diagonal de la matriz.
 */
int analizarTriangularizacion(struct MATRIZ *matriz) {

    // Recorremos la matriz por su diagonal
    int i;
    for (i=0 ; i < matriz->dimension ; i++) {

        // Si el elemento de la diagonal principal es 0, no es triangularizable
        if (matriz->elementos[i][i] == 0) {
            return 0;
        }    
    }

    // Ningun elemento de la diagonal es 0, por lo que es triangularizable
    return 1;
    
}



struct RESULTADO_FACTORIZACION triangularizar(struct MATRIZ *matriz) {
    
    // Creamos las matrices L y U.
    struct MATRIZ *matrizL = crearMatriz(matriz->dimension);
    struct MATRIZ *matrizU = crearMatriz(matriz->dimension);

    // Buffereamos los elementos de L y U, y su dimensión
    double **tempL = matrizL->elementos;
    double **tempU = matrizU->elementos;
    int n = matriz->dimension;

    // Inicializamos ambas matrices
    int i,j;
    for (i=0; i < matriz->dimension ; i++) {

        // Inicializamos U
        for (j=0; j < matriz->dimension ; j++) {
            tempU[i][j] = (matriz->elementos)[i][j];
        }

        // Inicializamos L
        tempL[i][i] = 1;
        for (j=i+1; j<n; j++) {
            tempL[i][j] = 0;
        }

    }

    double multiplicador; // Buffer del multiplicador por fila
    int k; // Iterador por elementos
    
    // Recorremos la matriz por columnas
    for (i=0; i<n; i++) {

        // Ahora debemos obtener 0 en la columna i, por debajo de la fila i.
        // Para ello usamos la fórmula: 0 = a - mx  <=> x = a/m donde:
        // -> x es el multiplicador q quiero encontrar
        // -> a es el elemento actual en la posicion de columna i
        // -> m es el elemento (i,i)
        
        for (j=i+1; j<n; j++) {
            
            // Calculo el multiplicador del elemento i,j.
            multiplicador = tempU[j][i] / tempU[i][i];

            // Modificamos la fila j
            for (k=i; k<n; k++){
                tempU[j][k] -= multiplicador * tempU[i][k];
            }

            // Guardamos el multiplicador
            tempL[j][i] = multiplicador;
            
        }
    }

    // Construimos el objeto a devolver
    struct RESULTADO_FACTORIZACION resultado;
    resultado.matrizL = matrizL;
    resultado.matrizU = matrizU;

    // Retornar el resultado
    return resultado;
    
}


/**
 * Este método toma como parámetro una matriz y encuentra su inversa asumiendo
 * que es una matriz triangular inferior.
 */
struct MATRIZ *invertirMatrizL(struct MATRIZ *matrizL) {

    // Buffereamos la dimension de la matriz
    int dimension = matrizL->dimension;

    int i,j; // Iteradores para obtener la inversa
    struct VECTOR *b = crearVector(dimension); // Buffer para la sustitición
    struct VECTOR *x = crearVector(dimension); // Buffer para la sustitición
    
    // Aquí guardaremos la inversa de la matriz
    struct MATRIZ *inversa = crearMatriz(dimension);

    // Obtenemos la inversa por filas
    for (i=0;i<dimension;i++) {
    
        // Armamos la columna b
        for (j=0;j<dimension;j++) {
            b->elementos[j] = (j==i) ? 1 : 0;
        }
    
        // Resolvemos el sistema
        int sustitucionExitosa = sustitucionDirecta(matrizL,b,x);

        // Si devolvió 0, no se pudo invertir
        if(!sustitucionExitosa) {
            return NULL;
        }
        
        // Copiar la solucion parcial a la matriz resultado
        for (j=0;j<dimension;j++) {
            inversa->elementos[j][i] = x->elementos[j];
        }

    }

    // Eliminamos los buffers
    eliminarVector(x);
    eliminarVector(b);

    // Devolvemos la inversa calculada
    return inversa;
    
}


/**
 * Este método toma como parámetro una matriz y encuentra su inversa asumiendo
 * que es una matriz triangular superior.
 */
struct MATRIZ *invertirMatrizU(struct MATRIZ *U) {

    // Buffereamos la dimension de la matriz
    int tam = U->dimension;

    int i,j; // Indices de iteración
    struct VECTOR *b = crearVector(tam); // Buffer para el cálculo de la inversa
    struct VECTOR *x = crearVector(tam); // Buffer para el cálculo de la inversa

    // Aquí guardaremos la inversa
    struct MATRIZ *inversa = crearMatriz(tam);

    // Calculamos la inversa por filas
    for(i=0;i<tam;i++) {
    
        // Armamos la columna b
        for (j=0;j<tam;j++) {
            b->elementos[j] = (j==i) ? 1 : 0;
        }
    
        // Resolvemos el sistema
        int sustitucionExitosa = sustitucionInversa(U,b,x);

        // Si devolvió 0, no se pudo invertir
        if(!sustitucionExitosa) {
            return NULL;
        }
        
        // Copiar la solucion parcial a la matriz resultado
        for(j=0;j<tam;j++) {
            inversa->elementos[j][i] = x->elementos[j];
        }

    }

    // Eliminamos los buffers
    eliminarVector(x);
    eliminarVector(b);

    // Devolvemos la inversa calculada
    return inversa;
    
}


/*
 * Este método realiza una sustitución directa para resolver el sistema
 * especificado en los parámetros.
 * Si no se puede realizar la sustitución, se devuelve 0.
 */
int sustitucionDirecta(struct MATRIZ *coefs, struct VECTOR *ti, struct VECTOR *solucion) {
    
    int i,j; // Variables de iteración
    double temp; // Buffer

    // Si el primer elemento de la matriz de coeficientes es 0,
    // no se puede determinar la inversa.
    if(coefs->elementos[0][0] == 0) {
        return 0;
    }

    // Obtenemos la primera incognita.
    solucion->elementos[0] = ti->elementos[0] / coefs->elementos[0][0];

    // Resolver las siguients incógnitas.
    for (i=1; i < coefs->dimension; i++) {

        // Si el coeficiente de la incógnita es 0, no se puede determinar la inversa
        if(coefs->elementos[i][i] == 0) {
           return 0;
        }

        temp = 0; // Reiniciamos el buffer

        // Calculamos el valor de la incógnita actual
        for (j=0; j<i; j++) {
            temp += coefs->elementos[i][j] * solucion->elementos[j];
        }

        // Asignamos la incónita calculada
        solucion->elementos[i] = (ti->elementos[i] - temp ) / coefs->elementos[i][i];

    }

    // Pudimos calcular todas las incógnitas
    return 1;
    
}


/*
 * Este método realiza una sustitución inversa para resolver el sistema
 * especificado en los parámetros.
 * Si no se puede realizar la sustitución, se devuelve 0.
 */
int sustitucionInversa(struct MATRIZ *coefs, struct VECTOR *ti, struct VECTOR *solucion) {

    // Buffereamos la dimensión de la matriz
    int dim = coefs->dimension;

    int i,j; // Variables de iteración
    double temp; // Variable buffer

    // Si el coeficiente de la incógnita es 0, no se puede determinar la inversa
    if (coefs->elementos[dim-1][dim-1] == 0) {
        return 0;
    }

    // Obtenemos la última incógnita
    solucion->elementos[dim-1] = ti->elementos[dim-1] / coefs->elementos[dim-1][dim-1];

    // Resolvemos las siguientes incógnitas
    for (i=dim-2; i>=0; i--) {

        // Si el coeficiente de la incógnita es 0, no se puede determinar la inversa
        if(coefs->elementos[i][i] == 0) {
           return 0;
        }

        temp = 0; // Reseteamos el buffer

        // Calculamos el valor de la incógnita actual
        for (j=i+1; j<dim; j++) {
            temp += coefs->elementos[i][j] * solucion->elementos[j];
        }

        // Asignamos el valor encontrado
        solucion->elementos[i] = (ti->elementos[i] - temp ) / coefs->elementos[i][i];

    }

    // Pudimos calcular todas las incógnitas
    return 1;
    
}

